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ABSTRACT 


A computer program is described for the calculation of the zeroes of the 

associated Legendre functions, P , and their derivatives, for the calculation 

nm 

of the extrema of P and also the Integral between pairs of successive zeroes. 

nm 

The program has been run for all n, m from (0,0) to (20,20) and selected 
cases beyond that for n up to 40 . Up to (20, 20) , the program (written In 
double precision) retains nearly full accuracy, and indications are that up to 
(40,40) there is still sufficient precision (4-5 decimal digits for a 54-blt man- 
tissa) for estimation of various bounds and errors Involved in geopotential model- 
ling, the purpose for which the program was written. 
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I. INTRODUCTION 

This report describes a computer program for the calculation of data on 
the associated Legendre functions of the first kind. These data are useful in the 
estimation of bounds for truncation error In the spherical harmonic expansion of 
the geopotential, and also for the estimation of bounds on the coefficients in such 
an expansion. The application of the results of this calculation to these estima- 
tion problems is discussed in References 1 and 2. The accuracy requirements 
for estimation purposes are not very stringent, a few significant digits should be 
adequate. The program can operate up to degree and order 100 ; this limitation 
Is imposed by the dimensioning of various arrays and would be easy to change. 

The program has been run from 0,0 through 20,20 and appears to have accu- 
racy of 8 or 9 significant digits for this range of degrees and orders. Runs 
for degrees 30 and 40 with order zero indicate that one can probably run It to 
40, 40 with an accuracy of four significant digits. The accuracy can probably be 
significantly increased by implementing one or another of the suggested modifica- 
tions to the subroutine for finding roots. > 

In constructing the program, two formulations for the associated Legendre 
functions were implemented. In one, z = cos 0 , where 0 ts the polar angle of 
spherical coordinates, is the Independent variable. In the other, x ** sin 0/2 is 
the independent variable. These two variables are related by 

z = l-2x (1.1) 

find the corresponding associated Legendre functions are given by 

va /2 2 

(1-z ) • polynomial of degree (n-m)/2 In z 

2 for n - m even (1.2) 

(1-z ) * z • polynomial of degree (n-m-l)/2 

in for n-m odd 
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( 1 . 2 ) 


P lx) = [x(l-x)j ' ■ polynomial of degree (n-m) ill x 

nrn 

From Eqs. (1. 2), it would at first appear that the calculation must accommodate 

three cases; actually there are six cases, since the extrema of are found 

from the zeroes of the derivative of P with respect to its independent vari- 

nm 

able, and the derivative must be handled in different ways for m = 0 and m>0. 
In addition, there are seven special cases that must be handled separately (e. g, , 
one of these is P » constant, for which there are no zeroes, extrema, or inter- 
val integrals. 


The "Interval Integrals*,’ mentioned above and in the title, are the integrals, 

between successive zeroes of P , with respect to its Independent variable, 

nm 

From Eq. (1. 1) 


dz = - 2dx 


x = 0 

corresponds to 2*1 

M- 

II 

corresponds to z = -1 


The final print-out of the full set of calculations lists the zeroes of P and its 

nm 

derivative in adjacent columns and in increasing order relative to the variable used. 
The associated extrema and interval integrals appear in the third and fourth columns. 
Because of the correspondence of the endpoints of the interval of definition of P^ 
indicated in Eq. the results read from top to bottom in the z-formulation 

correspond to those read from bottom to top in the x-formulation. The magnitudes 
of the zeroes are related by Eq. (1. 1). The extrema should be identical. The in- 
terval Integrals are related by a factor 2 which comes from Eq. (1. 3) (not -2 , 
since the minus sign Ls compensated by an interchange in limits of integration 
as one transforms from one formulation to the other). 


\ 
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In any particular run of the program, one formulation or the other is se- 
lected by an input switch. The two formulations were implemented because it 
seemed likely that they might well complement one another and, as we shall see, 
this is indeed the case. In addition, check-out of the program was greatly 
facilitated. 


Several output options are available through another input switch. The 
general flow of the main program is as follows: 


1. Input and initialization, including selection of the formulation 
to be used. 


2. Calculate and print the coefficients of the polynomial parts 


of P 
nm 

Option: 


and 


P’ 
nm 


Terminate the program at this point and 
go to more input at 1 


, 3, Calculate the zeroes of P and P* 

nm nm 

Option: Print these zeroes and go to 1 

Option: No print; bypass 4 and go to 5 

4. Calculate extrema of P by evaluation at the zeroes of P' . 

nm nm 

Option: Print zeroes and extrema and go to 1 


5. Calculate the interval integrals using the zeroes of P nm » 

Option (only if 4 is bypassed): 

Print zeroes and interval Integral and 
go to 1 

6. Print zeroes of P and P* , extrema of P , and 

nm nm nm 

interval integrals in tabular form. 

7. Go to 1 (with exit if no more data available). 


Listings of the main program and all the subroutines are provided in the 
appendices. The remaining sections of the report describe the steps listed above 
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in greater detail, with references to line and statement numbers appearing in the 
listings. 

Section n contains a list of the Input parameters and a discussion of their 
various functions. The branching involved in the six cases mentioned earlier is 
also described in Section H. This is followed, in Section HI, by the recursion 
formulas used to obtain the coefficients of the polynomial parts of P flm and , 
and a discussion of the subroutines in which they are implemented. 

The zeroes of the polynomial parts of P and P* are calculated by 

nni nm 

Graeffe’s root squaring method, implemented in subroutine GRAEFF. Some in- 
teresting problems were encountered, and these problems and their resolution are 
described in Section IV. This subroutine presently limits the accuracy of the pro- 
gram, and hence the size of degree and order to which it can be applied. The re- 
sults of a few test runs are presented, and several possibilities for improvement 
of the accuracy are discussed briefly. 

The extrema of P are found by direct substitution of the zeroes of P* 
nm um 

into P and this is accomplished by subroutines FUNCT and EVAL, which are 
nm 

straightforward and easily followed from the listing. The interval integrals are 
calculated by Gaussian quadrature in subroutine GAUSS, which is also straight- 
forward. A few comments on these three subroutines appear in Section V. 
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II. INPUT, INITIALIZATION, AND OUTPUT 

The major portion of the Main Program is taken up by Input, initialization, 
and output. The calculations are all done in subroutines, called by the Main Pro- 
gram. A listing of the Main Program is given in Appendix A. The references to 
symbols, statement numbers, and line numbers In this section apply to the Main 
Program. The output section is located between Statements 600 and 800. It fol- 
lows the flow indicated in the Introduction with the indicated options implemented 
in Lines 20800, 24600, 24800, 31300, and 31800. 

A block of 20 integers, IN(20), is reserved for Input parameters. A 
block of 100 integers, NUM(100), is also used for input under certain conditions. 
These blocks are in NAMELISTS INI and IN2 . The output of the program la 
carried In the arrays 


C(101) 

coefficients for the polynomial part of 

P 

nm 

CP(102) 

coefficients for the polynomial part of 

P T 

nm 

Z(102) 

zeroes of P 

run 


ZP(101) 

zeroes of P f 
nm 


EX(101) 

extrema of P 

nm 


FIN(lOl) 

Interval Integrals 



The first part of the initialization consists of identifying the input block, 
IN, with mnemonic names as follows: 

IN(1) = IND = 0 independent variable is z = cos 9 

2 

1 independent variable is x - sin 0/2 
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IN(2) = NOPT = 0 


EN(3) 

IN(4) 


IN(5) * 

IN(6) = 

IN(7) = 
IN( 8) = 
IN<9) = 

IN(10) 

IN(11) 

IN(12) 

IN(13) 

IN(14) 


a range of degrees equally spaced is desired; 
see IN(7), IN(8), and IN(9) 

> 0 a list of NOPT degrees to be read into the 

block NUM, using NAMELIST IN2 for the input 


MOPT = -1 
* 0 


process all orders consistent with each specified 
degree 

process only order MOPT for the specified 
degrees 


INC: Print Options: 

0 compute and print only C and CP 

1 compute and print only C , CP , Z , and ZP 

2 


! ITMAX 


NI 


compute and print only C , CP , Z , ZP , 
and FIN 

compute and print only C , CP , Z , ZP, 
and EX 

compute and print C , CP , Z , ZP , EX , 
and FIN 

maximum number of iterations allowed in 
GRAEFF for the calculation of Z and ZP 


use the zeroes and weight factors for P 


in GAUSS 


(NI+X) , 0 


IMIN 

: ISTEP > 
! I NX 

= NTOL 


IN(20) 


\ 


process a range of INX 
degrees starting at IMIN and 
spaced at ISTEP intervals 


convergence criterion 
SCALE = IN(1I)**IN(12) 

TOL = 10**IN(13) See Section V on GAUSS 

not used at present 


> See Section IV on 
GRAEFF 
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A single error return Is provided for several Input conditions which might result 
In poor functioning of the program. 

The second part of the initialization involves setting up the array NUM(l) 

th 

in such a way that NUM.(I) is the 1 degree to be processed, with a total of INX 
degrees. This information goes into the main DO loop starting at Statement 44; 

DO 1000, 1=1, INX followed by N1 = NUM(I), where N1 is the degree currently 
being processed. For NOPT>0, NUM is filled from the second READ statement 
(Line 4400). The DO loops to 6, 8, and 10 rearrange the degrees read and restore 
them to NUM so that 

NUM(Il) > NUM(I2) if and only if Il> 12 

This means that the degrees may be in any order in the data statement. For 
NOPT = 0 , Statements 20 and 30 construct NUM so that 

NUM(l) = IMIN 

NUM(I) - NUM(I-l) + ISTEP 

NUM(INX) « IMIN + ISTEP*(INX-1) 

Note that the dimensions of 101 and 102 for C and CP imply that the degree 
N1 must not exceed 100 . For direct input (NOPT>0) no test is made, but for 
NOPT = 0 , NUM(I) is not permitted to exceed 100 (see DO loop 30). 

The third part of the initialization calls subroutine FNORMO (Line 8600); 
this step, together with the call to FNORM in Statement 58, is better discussed in 
the next section dealing with the calculation of the coefficients in the polynomial 
parts of P and P ? . 

The fourth part of the initialization sets up IMX and the array MUM , 
which do for orders what INX and NUM do for degrees. If MOPT>0, MUM(1) = 
MOPT and IMX , the number of orders to be processed is set to 1 . If MOPT<0, 
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MUM and IMX are defined (inside the DO 1000 I = 1,INX loop) to include all 
orders consistent with the current value of N1 by the DO 45 loop. 


The final step in the initialization is perhaps the most complex; it starts 
at Line 10300 near the beginning of the DO 999 loop (which processes all orders 
specified for the current N1 value) and extends to Line 20400, just before CALL 
COEF. This step sets up the branching procedure for the six cases mentioned in 
the Introduction. A basic reason for the large number of cases was the desire to 
make use of the symmetry involved in the z = cos 0 formulation to reduce com- 
putation time. In this formulation, the polynomial parts of P and P' are 

nm nm 

polynomials in cos 0 , so that only their positive zeroes need be calculated and, 
from these, only the corresponding extrema and interval integrals need be cal- 
culated. The complete set ts then obtained from multiplication of this set by + or 
-1 . Further, there is little point in making GRAEFF find a zero root which is 
readily found by factoring. 


The parameter KIND identifies the six cases, the special case for each, 
and the differences in their treatment. The various parameters listed with KIND 
are as follows: 


NR 

NRP 

NC 

NCP 

NP 

NPP 


number of zeroes of P to be found by GRAEFF 

nm 

number of zeroes of P T to be found by GRAEFF 

nm 

number of coefficients in the polynomial part of P 

nm 

number of coefficients in the polynomial part of P T 

nm 

number of zeroes of ^ nrn > including endpoints and 
zero, if present 

number of zeroes of P 1 , including endpoints and 
zero, if present 


Parameters starting with K are used in the rearranging and augmentation pro 
cesses listed for each value of KIND below. 
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The case a = m = 0 is very special; there are no roots, extrema, or 
interval integrals. A special printout is provided as soon as this can be detected, 
Line number 9200. 

For m>L p T has zeroes at ±1 in the cos 0 formulation and at 0 
2 

and 1 in the sin 0/2 formulation; these points correspond to zeroes of P nm » 
rather than to extrema, at least for the purposes of this report. These zeroes of 
are ignored in the program and output. 

For IND - 0 (cos 0 formulation), most of the zeroes of P and P f 
' nm run 

are obtained by taking ± the square root of the output of GRAEFF. This formu- 
lation consists of four cases, as follows: 


KIND » 1: m = 0 , n even , special case is n = 2 

set of zeroes of P f must be augmented by ZP = 0 
nm 

extrema corresponding to zeroes of P* are sym- 
. nm 

metric about Z = 0 ~ 

v 

interval integrals are also symmetric about Z - 0 

set of interval integrals must be augmented by 
«lst zero 

X, and J 

x last zero 

KIND - 2 : m =* 0 , n odd , special case is n - 1 

set of zeroes of P must be augmented by Z = 0 
nm 

extrema and interval integrals are antisymmetric 
about Z - 0 

set of interval integrals must be augmented by end- 
point integrals 

KIND - 3: m^O, n-m even, special case is n = m 

set of zeroes for P must be augmented by Z » ±1 
nm 

set of zeroes for P' must be augmented by ZP = 0 
nm 

extrema and interval integrals are symmetric about Z = 0 
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KIND =4: m>0, n-m odd, special case is n=m+l 

set of zeroes for P must be augmented by Z = 0 , ±1 
run 

extrema and Interval integrals are antisymmetric about 
Z = 0 

For IND = 1 {sin** 0/2 formulation), subroutine GRAEFF gives all zeroes 

for P and P T except at the endpoints. The parity of n-m is not significant 
nm run 

and we do not exploit the symmetry properties of P^^ and P^ m about the point 


KIND -5: m = 0 , special case is n = 1 

output of GRAEFF is used unchanged for zeroes and 
extrema 

set of interval integrals must be augmented by the 
endpoint integrals 

KIND = 6: m>0 , special case is n-m 

set of zeroes of must be augmented by x = 0 , 1 

Although not properly a part of initialization, we mention here that in State- 
ments 140-220 the positive square roots of the output of GRAEFF are taken (for 
KIND = 1 ,2 ,3 ,4) and Z - 0 , ZP = 0 are introduced where necessary. The re- 
maining rearrangement of all roots, extrema, and interval integrals for output 
purposes is carried out In Statements 535-600. 

The special cases, identified by ISP - 1 , together with KIND , are given 
special treatment in Statements 800-910. 

A word should be said about values to be used for some of the input parameters. 

The principal reason for including ITMAX in GRAEFF was to avoid being trapped in a 

\ 

loop, in case convergence fails. The test cases run indicate that a reasonable value 
for ITMAX Is 
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ITMAX = 20 


since iterations in excess of 20 appear to have no significance. NTOL and 
SCALE are defined by IN(10), IN(U), and IN<12). The values used In testing 
the program were 


IN(XO) * NTOL = 14 


IN(ll) = 10 
IN (12) = 1 


implying SCALE = 10 


Utilization of a hexadecimal basis for SCALE with proper adjustment of NTOL 
might have computational advantages on the IBM 360. 


In all the tests carried out, we set 

NI - 9 

Some experimentation might show that a lower value could be used, particularly 
for small values of N, without sacrificing accuracy. Since there are NI+1 
evaluations of the integrand for each entry into GAUSS, some saving of machine 
time could be achLeved if lower values of NI yield acceptable results. In the 
tests on the program, we set 

IN{13) = -12 Implying TOL « 10 12 

This parameter is probably not significant for the analysis of P ; it was intro- 
duced so that GAUSS would be a self-contained subroutine, available for any program 
in which a Gaussian quadrature would be of use. 
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in 


CALCULATION OF THE COEFFICIENTS 


The coefficients for the polynomial parts of P and P' are calculated 

nm nm 

in three steps for the z = cos 0 formulation, using subroutines FNORMO, FNORM, 

2 

and COEF (listings given in Appendix B). For the x = sLn 0/2 formulation, only 
FNORM and COEF are required. We start by writing in the two formulations 

as 


with 


and 


P (z) 
nm' 7 


P (x) 
nm' 7 


m/2 C(n-m)/2] 

V-* 2 ) [I C nm (k > * ( “ >]: IND = 0 

k=0 


- (x(l-x)) 


m/2 


n-m 


C (k) x 
nm' 


k=0 


IND = 1 


(3.1) 


C (0) 
nm' 7 


= A 


nm 


C (0) = A 
nm' 7 nm 


C (k+1) = 
nm' 7 


C <k+l) - 
nm 


(n-m-2k)(n-m-2k-l) r "-m-2 

2(k+l)(2n-2k-l) urn 1 1 2 J 


1 T^£Tr t> 


nm’ 


(3.2) 


= jm! / ( 2 - 6 m o>< 2n+1 > 

nm 2 n *nl ' (n-m) ! (n+m) I 

A = —- 7 ( 2 - 6 1 )(2n+l) (3.3) 

nm ml V ' mO ' 7 (n-m) ! ' ' 

fl m0 = 0 m>0 • * [p] = larg6St integer S P 
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The C (k) appearing here are, of course, not geopotenttal coefficients. The 
nm , " r 

factors A and A include the factor 
nm nm 


i<2-S m o)( 2n+1) 


(n-m) ? 
(n+m) I 


( 3 . 4 ) 


which converts conventional associated Legendre functions into the fully normalized 
form used by geodesists. The derivatives for the two formulations of P^ m take 
the form, for m>0: 


dP 


nm 


dz 


/I 2 v 

(1-z ) 


<(m/Z)-l) 


[ 

LI 


n-m+l/| 


2 J 


CP (k) z 
nm' ' 


(n-m+l-2k) j 


k=0 


(3.5) 


dP 


nm 


dx 


((m/2)-l) “; m+1 
(x(l-x)) ) 

Li 

k=0 


CP (k) x 
nm' 


with 


CP (0) = B 
nm' ’ nm 


= -n A 


nm 


CP (0) = B 
nm nm 


m — 

-r* A 
2 nm 


CP (k) 
nm' ' 


C nm (k) - (n 2 -m 2 ) Cii _ 1>m (k-l)/(n(2n-l)) 

-m+11 
2 J 


= 12 T- 

A " t * * • * 1 


CP (k) 
nm 


t( m+ 2k)( t ‘ 2+n ) -m(m+k)(m+k-l)] 


k = 1 , 2 , . . . , (n-m+1) 


(3.6) 


Verification of these formulas is tedious, but straightforward. For m = 0, 
things are simpler: 
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\ 

\ 




CP n0 » . 


[n-m-l-2k] 


Hp n-m-1 

- 2-1 S nO«- k 


c JO) a B A » n A A 
nO' ' nO nO 


5 n 0<°> = B nO = A nO = ^ 


n —2k 

CP ft (k) = C A (k) 

nO' n nO' 


CP n0 (k) - (k+l)C n0 (k+l) 


The first step in the calculation of the C's and CP T s for the z = cos 0 

formulation is to calculate B A . This is done in subroutine FNORMG by setting 

no 


B oo - 0 


B 10 = S* 

and using the recursion relationship 


4k 2 - 1 

k-1 .k-1,0 


FNORMQ is called only once during a run, and computes and stores B up to 
and including the maximum value of n to be processed. Is so simple that 

it is calculated when needed in subroutine FNORM. 


The second step In calculating the coefficients is carried out in subroutine 

FNORM, which computes A , B or A m , B^ f depending upon the formu- 

nm nm ntn nm 

lation selected. For m = 0, of course, the calculation is trivial. For m>0, 
the following recursion formulas are implemented in FNORM. 


nm 


hi 


- / 

I 2n 

V n+1 


n-m+l 

— ; — A , 
n+m n, m-1 


— A 
n+1 nO 


m> 1 


A n0 = V" 


B 


nm 


- - n A 


nm 


n,m+l 


_ (n+m+l)(n-m) r* 


m+1 


nm 


m>l 


(3. 11) 


nl 


= J 


2{2n+l)(n +n) 


A _ = 2n+l 
nO 


B 


nm 


* A 
2 nm 


The factor (2-6 J in' A and A necessitates starting the recursion from 
ml) nm nm 

A_, and A^, , rather than from A^ and A nn , 


nl 


hi 


nO 


Finally, subroutine COEF, using the output of FNORM, implements the 
recursion formulas given in Eqs. (3.2), (3.6), and (3. 8) to obtain the C's and 
CP's or C's and CP's , depending upon the formulation desired. 


No study of the growth of error with the number of passes through these 
recursion formulas has been made. It has been noted by S. Pines (Ref. 3) that 
care must be exercised in the use of recursion formulas. It is possible that in- 
accuracies in the coefficients are responsible for the lack of precision in the 
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determination of the zeroes of P and P f f although the way In which this 

nm nm 

occurs suggests that other effects dominate any inaccuracy in the coefficients. 
This matter is discussed further in the next section. 

Note that slight variations appear between the formulas given in this sec 
tion and their Implementation in subroutines FNORMO , FNORM, and COEF, be 
cause DO loops cannot start from zero. 
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IV. THE GRAEFFE HOOT SQUARING METHOD 


The subroutine for finding the zeroes of P and P 1 Is GRAEFF (AA f 

am nm 

N, Z, SCALE, NTOL, ITMAX, IND). The lLstlng is in Appendix C. It calculates 
the zeroes of a polynomial of degree N-l , with a coefficient array AA of M 
elements, associated with increasing or decreasing powers of the variable according 
as IND .is 1 or 0 . The Graeffe root squaring method is implemented in less than 
full generality: An implicit assumption is that the roots are real, positive, and dis- 
tinct, a condition fulfilled by the polynomial parts of P and P f , if z Is fac- 
tored from those of odd degree In the cos 8 formulation. The zeroes are stored 
In the array Z. The remaining entries in the calling sequence, SCALE, NTOL, 
and ITMAX will be discussed later. 


First we outline the basic idea of the method; an excellent dLscussion is 
given by Lane zos (Ref. 4), We suppose that 


x;>x > . . . >x >0 
12 n 


(4. 1) 


are the zeroes in descending order of magnitude of the polynomial 

n 


k=0 


(4.2) 


Then 


2 2 2 
y l “ X l >y 2 = *2 > " ■ • >y n = *n 


(4.3) 


are the zeroes in descending order of magnitude of the polynomial 

n 


I v‘ 

1=0 


(4.4) 
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where 


B. = a: 


B. 


A 1 + 2A 0 A 2 


B. 


A 2 ” 2A 1 A 3 + 2A 0 A 4 


B 


n-1 


(-l) n_1 (A 2 - 2 A A ) 

' ' ' n-1 n-2 n' 


* n 2 
(- 1 ) A n 


(4.5) 


As this process Is Iterated one obtains, on the K iterate, a polynomial with 

[k] 

coefficients B and zeroes 


< 2K )>*< 2K >> 



(4.6) 


th st 

such that the ratio of the i to the (i-1) zero becomes arbitrarily small 

for all i and sufficiently large K . Using this fact, and the relationship between 

[ K ] 

the coefficients B and sums of products of roots, it is easy to verify that 



(4.7) 


(or its reciprocal, depending on IND) converges to the zeroes of the given poly- 
nomial, As the iterates of the coefficients are constructed, It becomes ap- 
parent that they become more and more widely separated in order of magnitude. 
Numerically, the method terminates when the separation of the coefficients be- 
comes such that 





(4.8) 


IB 



because the remaining cross-product terms [see Eq. (4,5)] are beyond the word 
length of • If the word length for the calculation is L decimal digits, 

the criterion for termination is thus 

2B [ + f- B H 3 <[ B S IC3 ] 2 - 10 ' L 

for all relevant values of j , This is essentially the criterion used in GRAEFF, 
and L is given the name NTOL , an input quantity. 


In this subroutine, the terms contributing to each are added on one 

at a time from left to right, as shown in Eq. (4, 5), An array K1(I) is defined 

[ K ] 

to give the number of terms making up B. from the previous set of coefficients 

[K-l] 1 

. When the last term in this sum is beyond the word length of the K1(I)-1 

terms already summed, K1(I) is diminished by 1 . When K1(I) - 0 for all I , 

the iteration terminates. 


Since both round-off error and machine time can be expected to increase 
with the number of iterations, ITMAX , another input quantity, is also allowed 
to terminate the iteration, in which case the calculation of the zeroes proceeds 
on the basis of the B's so far obtained. In this case, a message is written to- 
gether with the array K1(I) , which indicates which of the B’s have failed to 
converge. An error message is written if any zero is negative, and the calcula- 
tion proceeds with the absolute value of such a zero. A standard print states the 
number of iterations used on the current entry to the subroutine. 

A significant problem in the implementation of Graeffe f s method arose be- 
cause the iterates of the coefficients grow very rapidly, and soon produce overflows. 
To avoid this problem, the parameter SCALE is used to convert all coefficients and 
their iterates to values less than SCALE and greater than or equal to 1 . Then 
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additional arrays are introduced to carry the powers associated with the coefficients; 
i« e, , for each I 

1 * B(I) < SCALE NEXB(I) = power (4. 10) 

and the actual corresponding coefficient is given by 

B{I) * SCALE ** NEXB(I) (4.11) 

The program has been run (in double precision) using SCALE = 10 , 

NTOL - 14 for all orders and degrees of P from 0,0 to 20,20, on the 

nra 

DEC KA10 , which has a mantissa of 54 bits. The indications are that the 

zeroes near zero hold 15 decimal digit precision for polynomials at least up to 

degree 20. The polynomial parts of and have their largest zeroes 

near unity and for such a polynomial part of degree 10 , the largest zeroes have 

10-11 digit precision; for one of degree 20 , the precision of the largest zeroes 

is only three or four digits. These data on precision were obtained by comparison 

of the zeroes of P ^ tabulated by the National Bureau of Standards (Ref. 5 ) , and 
nO 

by comparison of the output from the two formulations. In fact, the availability 

of the two formulations probably enables one to go to 40 ,40 with 7-8 digits of 

2 

precision. This is so because the small zeroes of the sin 0/2 formulation can 
be transformed into the zeroes near unity of the cos 0 formulation, while the 
small zeroes of the cos 0 formulation are transformed into those near x = % 
in the sin 0/2 formulation. Thus, using the ’’good" zeroes from each of the 
two formulations and the symmetry properties, a set of zeroes good to 12 or 13 
digits may easily be constructed for 20 a 20 . Selected cases up to 40 , 40 have 
been run. The user of the program is cautioned that an ITMAX of 20 will be 
exceeded and that overflows may occur for IND = 1 , if n-m is appreciably 
greater than 20. Both conditions may be ignored since they affect only those 
zeroes for which significance is already lost; they must be found by running 
IND = 0. 
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One would like, of course, to account for the lack of precision of the "large” 

zeroes and, if possible, Improve the accuracy. An immediate thought might be 

that errors in the input coefficients (recall that these are computed by recursion 

formulas) are the primary cause. This does not seem likely, however, because 

for IND = 0 (cos 0 formulation) the most important coefficients for large zeroes 

are those which start the recursion calculation. Still, all coefficients do ultimately 

fKl 

enter the iterates of for the large zeroes, and the possibility cannot be 

eliminated without further testing. Another thought is that the round-off error 
produced by scaling Is the culprit. This possibility has been tested and round-off, 
while present, is several orders of magnitude less than the discrepancies observed. 
The most plausible, but as yet untested, explanation is loss of significance in the 
subtractions implied by Eq. (4. 5). It is possible that combining these terms starting 
with the smallest and ending with the largest (in magnitude) might help, at the ex- 
pense of machine time spent in the sort. Probably the most practical method to Im- 
prove the situation is to use the output of GRAEFF as the initial guess to a Newton 
procedure, ^ 
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V. 


CALCULATION OF THE EXTREMA AND INTERVAL INTEGRALS 


The subroutines for these calculations are straightforward and require little 

comment. To obtain the extrema, P must be evaluated at the zeroes of P' . 

nm mm 

This calculation is carried out by subroutines FUNCT (X,F ,L) and EVAL (A ,N, 

M , X , P , IND) (listed in Appendix D). Subroutine EVAL simply evaluates a poly- 
nomial of degree N-l with a coefficient array A (associated with ascending or 
descending powers of the variable, according as IND = 1 or 0) at an array X of 
M points. These evaluations are returned in the array P. Subroutine FUNCT, 

which accepts the array of L evaluation points X , supplies whichever of the 
o m/2 2 m / 2 r -m/2 

factors (1-z ) , z(l-z ) , or Lx(l-x)j is applicable and returns the 

values of P in the array F . It appears that the extrema are relatively Lnsen- 
nm 

sitive to errors in the zeroes of P* . This supports the opinion given in the pre- 

nm 

vious section that errors in the coefficients C and CP are relatively unimportant. 

To obtain the interval integrals, subroutine GAUSS (A , B ,NI , ABINT ,TOL) 
implements the Gaussian quadrature procedure, which is well described by Lanczos 
(Ref. 4). Two input options are provided: The zeroes and weight factors for P^q, 
k = 2 , 3 , . . . , 10 , are stored in data statements. The parameter NX selects 
those for k = NI+1 . A parameter TOL is introduced to avoid difficulties with 
small differences: If the limits A , B of the integral to be evaluated satisfy 

|a-b|<tol (S.i) 

the subroutine returns the value zero in the output parameter ABINT , and prints 
out a message to this effect. Subroutine GAUSS calls FUNCT and then EVAL to 
evaluate the integrand where necessary. 



The interval integrals are quite sensitive to inaccuracies in the zeroes of 

P , as one might expect, since these Inaccuracies will destroy the non-negative 
nm 

character of the integrand. However, it is felt that, using both formulations and 
symmetry considerations, the interval integrals have 8-10 digits of accuracy up 
to 20 , 20 and will probably retain 3-4 digits perhaps up to 40 , 40 , which 
should be adequate for the estimation purposes discussed in Reference 2. 

It should be mentioned that the program does not implement the construction 
of a single table for the zeroes, extrema, and interval integrals utilizing the output 
of the two formulations in such a way as to maximize accuracy. The necessary 
additions to the program would be easy to Insert, Time, however, did not permit 
sufficiently detailed examination of the output to determine the points at which the 
switch between formulations should be made. These switch points are very likely 
functions of m and n f though perhaps sensitive only to the difference n-m. 
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WRITE (6,620) UP(K) jKsl,NPP) 

FORMAT (* HEROES OF PNM PRIME APE 1 f / < 10X , 1P3D23 , J4 ) / > 

IF (INC * 2) 999,650^630 
WRITE (6,643) (LX(K) ,K31,NPP> 

FORMAT {! EXTREMA OF .PNM ARE \ ' /< J.0X, XP.3D25, 3,4 )/) 

IF ( I N C E 0 , 3 ) GO TO 999 
WRITE (6,663) (FtN(K)iKsliNFP) 

FORMAT (# INTERVAL INTEGRALS FOR PNM ARE I * f 
C(li:X, iP3D25 t l 4 )/) 

GO TO 999 
WRITE (6,713) 

FORMAT (* HEROES or PNM ANO PNM PRIME, 
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WRITE (6,730) HP(K),EX(K),F|N<KK) 

FORMAT (20X,XP3O25 | X4) 

K2eK+i 

WRITE (6,740 ) H ( K 2 ) 

FORMAT (3X, 1PU25,14) 

CONTINUE 

IF (M1.GT.0) GO TO 999 
WRITE (6,760) FJN(NPi-l> 

FORMAT (43X,f INTEGRAL FROM LAST HERO TO X IS 
Cl * ,4X,1PD25 I 14> 

GO TO 999 

WRITE (6,770) 

FORMAT (l P 00=1,0) P00 PRIME *0 | NO ROOTS, 

AMO EXTREMA, NO INTERVAL INTEGRALS * ) 

GO TO 1030 

GO TO (810,620, 530,840,850,860) KJNO 
H(2) s DSQRT (*»C(2)/C(1) ) 

Z(l)s-2(2) 

HP(1)»2,O0 
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36 100 


cm)BC(2> 

36200 


CALL CAgSS<E<2> i 1 ,D0#NI 4 rjN(3) iTOL) 

363 fcjfe 


CALL GAUSS (E (1) |H(2? f N!,F IN<2 >iT0U> 

36400 


ri^U)*F!N<3) 

36500 


GO TO 603 

36600 

620 

H(l>«0 t O3 

36700 


CALL GAUSS(0', D3| 4,D0i NJ ,rjN(2> «TOL> 

36800 


FlV(i)»-nN(2) 

36900 


GO TO 903 

37P300 

930 

2(1)5*1,00 

373.00 


2<2>*liD3 

37200 


2P(l)s0,D0 

37300 


E!X ( 1 ) s C ( 1 ) 

37400 


CALL GAUSS{®l t D0|l,D£fl,NJ # FlNj[l),TOL) 

37500 


GO TO 603 

37600 

840 

2(1)3*1,00 

3 7 700 


2 ( 2 ) =0 , 03 

37800 


2 ( 3 ) s 1 , 03 

37900 


ZPC2) =0SQRT(cCP(2)/gP(l) ) 

38000 


?P(l)a-HPC2) 

38100 


CALL FUNCT(HP,EX,2) 

38200 


CALL GAUSS(3 I O0,1 I O0»NI,FIN(2),TOU 

38300 


FPM1) 5 -FJN(2) 

33400 t 


GO TO 603 

38500 

850 

2 ( 1 ) s *C ( 1 ) /C ( 2 > 

38600 


CALL GAUSS(3 1 O3,2(l)#NI l FlN(i)|T0L) 

3 8 7 0 0 


FIM(2)=rFJM(l) 

38800 


GO TO 900 

38900 

960 

2(1)50,00 

39000 


2<2)=1,03 

39100 


2P(l)s 1/2,00 

39200 


CX(X)sC(U/2',U0**Hl 

39300 


CALL GAUSS(0’,O3,1,O0#N1 ,FIN(i) # TOL) 

39400 


GO TO 603 

39500 

900 

IF ( K I N'D L Q , 2 ) N X 1 s * 1 

39 600 


IF (KlNO',tO,5)NXl»0 

3 9 700 


WRITE (6,910)Z(1) ,NX1,FIN(1) ,FIN(2) 

39800 

910 

FORMAT (f P10 HAb ONF 2ER0 AT *, 1PD13 , 2 , AND 

39900 


A EXTREMA', THE INTERVAL INTEGRALS ARE I ' /3*,' FROM* 

40000 


B # 13# ^ TO 2(1)1* #iPD2*,14,' M,5X,'FR0M 2(1) 10 

40100 


Cll ' ilP025,14) 

40200 


CO TO 999 

40300 

999 

CONTINUE 

40400 

1003 

continue 

40500 


GO TO 5 

4 06 0*0 

2000 

WRITE < 6 1 2010 ) 

40700 

2013 

FORMAT (' INPUT JN1 DEFECTIVE! CO TO NEXT CASE'} 

40800 


GO TO 5 

40900 

3003 

CONTINUE 

41000 


END 
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APPENDIX B - Listing of the Coefficient Subroutines 


00100 


SUBROUTINE TNO«H0 ( N# J NO) 

00200 


IMPLICIT REAL*® 

00300 


COMMON /NOHM0/BU0U 

00400 


IF UNO **1)101 100* 200 

00500 

10 

911)50,000 

00600 


B(2)=DSQRT (3,000) 

00700 


00 20 U2fN 

00800 


El«i 

00900 


EISQ=EI*EI 

01000 


ENUM = DSQRT (4’ f 0D0*EISU*1 ,000} 

01100 


8( I-M)HNUM*BU)/(U*1,O0) 

01200 

20 

CONTINUE- 

01300 

01400 

100 

RETURN 

RETURN i 

return V 

01500 

200 

01600 


END 





SUBROUTINE FNQRMIN.M* INQ) 

00200 


IMPLICIT REAL *8 {A»H,OpZ> 

00300 


COMMON /NORM0/ B0<101) 

00 4 00 


COMMON/NORM/AUW »BU8J1) 

00500 


EN *N 

00600 


IE < I NO **5- ) 10,100|Z00 

00700 

10 

A(1)=B0(N+1)/EN 

00800 


BU)sB0(N+1) 

00900 


IE ( M , EQ', 0 ) RETURN 

0X000 


I E ( M L T , 0 ) NM*N 

0X100 


1F(M', CT, 0) NMsM 

01200 


C s Z,O0*EN/(EN + 3,,O0 J 

01300 


IF ( M , EQ'i 1 ) CO TO 90 

01400 


3(2)=B(1)*BU>»C 

01500 


DO 30 I » Z|NM 

01600 


El = I 

01700 


D = (ENbEI+I'iOBI/IEN+EI) 

01800 


B(I+1)«B(D*0 

01900 

30 

CONTINUE 

02000 


IFtM’, GT,0> CO TO 00 

021C0 


AU)=BU>/EN 

02200 


NN*N+1 

02300 


DO 52 I s 3 1 NN 

02400 


0 ( ! ) s-DSORT (0U)| 

02500 


ACI>a«Bf I>/EN 

02600 

50 

CONTINUE 

■02700 


CO TO 90 

02800 

60 

BtM + l)CA»DSQRT<B<M+l) ) 

02900 


A(M + l)=*8(MU)/EN 

03000 


RETURN 

03100 

90 

C * OSQRT(C) 

03200 



03300 


A ( 2 ) ( 2 >/EN 

03400 


RETURN 

03500 

100 

EN2PN » EN*(EN"**1|O0) 

03600 


C * ( 2 » Q0*EN + 1 i 00 ) 

03700 

U0 

A ( 1 ) sQSQR T ( C ) 

03820 


B(1)s-EN2PN*AU) 

03900 


IF <M,EQ t 0) RETURN 

04000 

120 

A<2)=C*EN2PN«2,U0 

04100 


IF tM,LT|0> MM* N 

04200 


IF ( M", CT , 0 ) MHoM 

04300 


DO 130 \ * 2 i MM 

04400 


El = I 

04500 


E12MI * EI*EI*£l 

04600 


AUH) s AU )«(EN2PN^EJ2MJ >/(EI*en 

04700 

130 

CONTINUE 

04800 

150 

EM a M 

04900 


IF (MMiEQiM) GO TO 160 

05000 


DO 160 J *1 # MH 

05100 


El * I 

05200 


A ( I+l) a QSQRT(A< 1*1) ) 

05300 


R(1H) 8 { AU*in«EI/2,O0 

05400 

160 

CONTINUE 

05500 


RETURN 

05600 

180 

A(M + 1)**DSQRT ( A U1*1 ) ) 

05700 


B^n) a A{Mn)*£M/2,p0 

05800 

200 

RETURN 

05900 


END 
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00100 
00200 
00300 
00400 
00500 
00600 
00700 
00800 
00V00 
01000 
01100 
01200 
01300 
01400 
01500 
01600 
01700 
0 1 80 0 
01900 
02000 
02100 
02200 
02300 
02400 
02500 
02600 
02700 
02800 
02900 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
03900 
04000 
04100 
04200 
04300 
04400 
04500 
04600 
04700 
04800 
04900 
05000 
05100 
05200 
05300 
05400 
05500 
05600 
05700 
05800 
05900 
06000 


SUBROUTINE COEF 
IMPLICIT HEAI*«<A-H|0-H> 

C COMPUTE COEFFICIENTS OF POLYNOMIAL PARTS OF PNM # PNM PRIME 
C INPUT: ORDER N, DEGREE M AND 1ND TO GIVE FORMULATION DESIRED 
C OUTPUT \ NUMBER NC~OF COEFFICIENTS C OF PNM 
C NUMBER NCP OF COEFFICIENTS CP OF PNM PRIME 

c coeficients cm, cpu) 

c IND?0t PNM GIVEN IN POWERS OF ( COS ( THETA > ) #»2 WITH 
C NC*C (N«M+2>/2J# NCP= NC+i 

C INQ = i: PNM GIVEN IN POWERS OF ($ l N{ THETA/2 ) )t*2 WIIH 
C UC*N«M# NCPsN*M+i 

C IF Hs0 # nCPsNC« 1’F0R BOTH FORMULATIONS 
C JND'.GT.l MAY BE USED FOR OTHER FORMULATIONS 
CO M MON/NORM/A (101) ,8f 101) 

COMMON/COEFF/C(101) ,5P<102>|NC,NCPfN|M, INO 

NtfM=N*M 

enmm=nmh 

ENsN 

RMsM 

na)sA(M + l) 

CPU>sB<M + l) 

EN2PNsEN«EN 4EN 
IF( INO-l)100i?00 , 700 
100 ENMO2=ENMM/2,D0 

NMM0?sENM 02 
K 1 s NC« 1 

TWONPl«2',D0*eNn’,O0 
TWONM 1-TWO Npl**2'| 00 
fwOH=2,D0«EM 
ENMPlsENMM+l|D0 
T=-ENhO2»CP(1)/(EN»JW0NM1) 

3 = EN?PN«»EtN<*EM+TWQM 
DO 150 K»1 #K1 
EK = K 

TW0Ks2,D3«EK 

C(K+i ) = *C(K)*UENMM + 2, O0»TWOK)«<ENMP1pTWOK) 
1/(TW0KMTW0NP1*! w QK) > 

CP ( K«4*i } sT*S 

Ts-T*(ENMM«TW 0K)MENMP1*TW0K>/( (TWOK 
l + 2 t O0)MTWONMl"TWQK) ) 

SsS+TWOM 

150 CONTINUE 

IF CM’, EQ,0) CO TO 200 
I F ( NC , GE', NCP ) RETURN 
ENCaNC 

CP(NCP)=C(NC) 

RETURN 

200 CO 210 J 42i NC 

e I a I 

CP<I)s<ENMM*2,D0*(El*l,D0))*C(n 
210 CONTINUE 

RETURN 

400 UO 420 KeliNMM 

EKeK 

TW0K=2,D3*EK 

EMPK«EM*EK 

EMPKM1=EMPK«1,O0 

EKMPKsEK«EMPK 

T3(EN2PN"EMPKM1«EMPK)/EKMPK 
C ( K+l ) * w T#C (K ) 



06*00 


SstEM+TW0K>*EN2PN-EM*EHpKfU*EM p K 

06200 


CP<Kn>="5*CU>/<2fO#*EKMPK> 

06300 

420 

CONTINUE 

06400 


CPINCPJs-ENaCiNC) 

06500 


IF < M GT i 3 ) RETURN 

06600 


DO 450 K*2 # NCP 

06700 


ER C K 

06800 


CP ( K ) sEK«C < K + l ) 

06900 

450 

CONTINUE 

07000 


return 

07100 

703 

RETURN 

07200 


END 
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00100 
00208 
00300 
00400 
00500 
00600 
00700 
00800 
00900 
01000 
01*00 
01200 
01300. 
01400 
01500 
01600 
01700 
01800 
01903 
02000 
02100 
02200 
02300 
02400 
02500 
02600 
02700 
028 00 
0g900 
03000 
03100 
03200 
03300 
03400 
03500 
03600 
03700 
03800 
03900 
04003 
04100 
04200 
04300 
04 400 
04500 
04600 
04700 
04800 
04900 
05000 
05103 
05200 
05303 
05400 
05500 
05600 
05700 
05800 
05900 
06000 


APPENDIX C - Listing of Subroutine GRAEFF 

SUBROUTINE GRAETF ( AA,N,2#5CAlE,NT0i_# jTMAXi I ND > 

IMPLICIT REAl*8 ( A-H#0-H> 

DIMENSION A ( 102 ) , AA ( 102 ) i 2 ( 102 ) , 8 ( 102 ) #K 1 < 102 > I 
CNEXA ( 102 ) # N£XB < 1 ^ 2 ) 

C ROOTS 2(1) OF A POLYNOMIAL OF DEGREE N^i 0Y GRaEFFE‘S METHOD 
C Mi) IS COEFFICIENT OF X**(Nd)i lal«2i.«.N FOR I NO * fl 

c w « M of x**u-i). w for iNb.cT .0 

C ITMAX is THE MAXIMUM NUMBER OF ITERATIONS 
C KTOL TFRMINATES ITERATION ON A CONVERGENCE CRITERION' 

IF ( N * E Q * i ) GO TO 270 
IF ( N , E Q V 2 ) GO TO 28 0 
ITER ■ 1 
EN = N 

FN02 s EN.2,000 \ 

N02=EN02 

no 10 I 5 1 # n 

A(I)aAA(l) 

S1G1^1,0D0 

IF (i.LEV N 02) KlC !>■!»! 

If (i . gt,no2) ki < I ) =n» i 

NEXS0 

T E STs A < | ) 

IF (fESTVLT.0,00) SlGls-1,000 
TESTsOABS(TeST) 

2 IF (TEST.LT.SCALE) GO TO 4 

TEST 3 TEST/ScA L E 
NEX= NEX^l 
CO TO 2 

4 IF (TEST, GE. 1.00) GO TO 6 

tEST c TEST#SCALE * 

NEXsnEX * 1 > 

GO Th 4 

6 NEXA(I)=NEX 

A ( I ) a SI Gl*TEST 
le continue 
2 2 no 100 I C 1#N 

SIG=-1.0O0 

CsA(h«A(!> 

NEXCsNEXA(I)»2 

KsuM^KKI) 

IF (K1(I) • EO , 0 ) KSUM si 
3g DO 95 K s i , KSUM 

IF ( k 1 < I > , E q , 0 ) G 0 Tn 75 
fERM*A(I*K)*A(I-K)*2 ( 0D0 
NEXT,NEXA( UK)*NEXA< UK) 
nEXO* N EXC-nEXT 

IF( UBS(NEXO) .LTVNTOL) go TO 45 
IF (K.EQVkSUM) Kl ( I) »KSUM^i 
CO TO 94 

45 IF <nEXDVlT,0> Go To 50 

f ERM s TERM *SC ALE (»NEXD) 

C = C4.TERM«SJG 
GO TO 75 

52 C=C*SCALE#*(NEXO) 

C=C + TERM*SIG 
NEXCaNCXT 
75 S I Gl s l 

IF(CVlT,0 # O0)SIG1=«1 

C B OAB$iC) 

80 IF (CVlT.SCALE) G 0 To 85 
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06100 


C=C/*CALE 

06200 


NEXC'a NEXC*i 

063^ 


GO TO ©0 

06400 

83 

IF <C'aGE«l«O0) GO to 90 

06 500 


C-C*SCALE 

06600 


MEXCsnEX C*l 

06700 


GO TO ©5 

06800 

90 

C=C*S!Gl 

06900 


!F <KlCl>,£O'»0) GO TO 96 

07000 

94 

S iG a »$lG 

07100 

95 

CONTINUE 

07200 

96 

B 1 1 > B C 

073 0a 


NEXB ( I ) sNEXC 

07400 

180 

continue 

07500 


no 110 1 8 l » N 

07600 


IF (kKU.GrVO) GO TO 120 

0770fl 

llfl 

CONTINUE 

07800 


Go T 0 200 

07900 

120 

IF ( ITER i GE , 1 TMAX) GO TO 180 

08000 


I TER = I TER+i 

08i00 


8IG=i*OO0 

08200 


00 130 K=ltN 

08300 


a<k>=q<k>*sig 

08400 


NEXA (K)sNeXb<K> 

08500 


SlG=7SlG 

0 Q 6 0 0 

120 

continue 

08700 


GO TO 20 

08800 

180 

WRITE <6»l90> <Ki<I)Il=liN) 

089 0 fl 

19 0 

FORMAT r ITMAX EXCEEDEDI K ISM 0 l5/<3X#15l5/>> 

09000 

200 

rxpa?.n0**<-iTER) • 

09100 


n3»N-1 

09200 


00 250 I *1 # N3 

09300 


IF n'NO) 210# 210 * 220 

0 94 00 

2l 0 

n4 = n;i 

09500 


?1 

09600 


NEXH = NEXB ( I + 1)^NEXB( I ) 

09700 


GO TO ?3 0 

09800 

220 

?1 = fl< I )/BU*l> 

09900 


NEXH«NEXB( I ) ff NEXD< !♦£) 

10000 


N4 = I 

10100 

230 

IF <?1 i LT ♦ O ) WRITE < 6 1 240 > n4,Z1 

ifj20a 

240 

format r 2(M5 #») negative AND EQUAL T o f ; E i 8 V 8 J 

10300 


21 sDABS(Hi )#*EXp 

10400 


EX2=NEX2 

10500 


rx? = FX5?*ExP 

10600 


7(N4>=21 *ScALE**EX2 

10700 

. 250 

continue 

1 0. 8 0 0 


WRITE (6i260> ITER 

10900 

260 

FORMAT! * GRAEFF USED ♦ , 1 3 ; ♦ I TER AT I QN$ » ) 

11000 


RETURN 

11100 

270 

WRITE ( 6 # 2 7 5 ) \ 

11200 

275 

FORMAT <» POLYNOMIAL IS OF DEGREE 0 | NO ROOTS*) 

11300 


RET^N ... 

11400 

260 

IF UND.EQ.0) Z<i)*-AA(2i/AA<l> 

11500 


IF < I NO » Eq * l ) 2(l)#i!'AAU)/AA(2) 

11600 


W R j T E <6i205) 2(i> 

11700 

285 

FORMAT <» POLYNOMIAL IS LINEAR] 2(1) «»#lP025 l 14) 

11800 


RETURN 

11900 


END 
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APPENDIX D - Listing of GAUSS and the Function Evaluation Subroutines 

<40100 SUBROUTINE TuNCt 

00200 IMPLICIT REAL*fl<A»HfO*H) 

00300 C TO EVALUATE ASSOCIATED LECTURE FUNCTIONS PNM AT 

00400 C L POINTS X, OUTPUT JS L VALUES OF PNM, INO 

00500 C INDICATES THE FORMULATION USED*, 

00 600 COMMQN/COEFF/CUGU) # CP (102) ,NC,NCP,N|M, INO 

00700 DIMENSION XUtTO«Pu62)|Y(3,03) 

00800 EM = M 

00900 EM02=EM/2 ,0D0 

01000 IF(lND-l) 10 | 100 i 200 

01100 10 DO 20 1=1, L 

01200 V(I)sX(n«X(I) 

01300 22 CONTINUE 

01400 CALL E VAL ( C # NC # L i Y | P I NO > 

01500 IF (M0O<<N*M) ,2) ,£0,0) CO TO 40 

01600 RO 30 I =1 , L 

0172 0 PCI )sp< l)«X(l )*UiO0»X< J)*X( {))**EMQ2 

01800 30 CONTINUE 

01900 RETURN 

02000 40 DO 50 1 = 1, L 

02100 P(l >-sP( ! )»(1|QU*X( I )*X( J ) )**EH02 

02200 50 CONTINUE 

02300 RETURN 

02400 102 CALL EVAL<C,NC,L#X,P# JND) 

02500 00 lie i =1, L 

0 2600 P(I)3p(IJ»(XU)*(l,0tJ-XU) ) )««CHQ2 

02700 110 CONTINUE 

0H80P RETURN 

02900 200 RETURN 

03000 END 
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»* 3 I 1 '('■ V* 

Vy&\Z 

03300 

00400 

0050? 

03600 

00700 
a ^ b 0 55 

»rJ9fej{? 

01003 
01100 
012055 
01303 
01400 
01500 
01600 
0^730 
31800 
01 9GB 
02000 
02103 
02200 
02303 
02403 
02500 
026 00 
02700 
02800 
02903 
03000 
0 3100 

03?^ 

03300 

03400 


SUBROUTINE £VAL<A,N,M,X#P# I ND > 

C EVAIUATC A POLYNOMIAL OF DEGREE N«*l AT M POINTS XU) 

C return M values PCI) 

C A ( I ) IS THE COEFFICIENT OF X#«(M-I) FOR JNQ * 0 
C « " »» « *! 5 1 1 FOR I NO » NE » 0 

IMPLICIT RE AL*8 (A-H,0-2> 

DIMENSION AU02),X(102)#P<102> 

IF < N , GT , 1 ) GO TO 10 
DO 5 K«1>M 
P(K>=A<1> 

5 CONTINUE 

return 

10 IF < I NO | EQ , i ) CO TO 52 

no 40 1 = 1 , M 

T=X(I) 

Y = A< i ) *T*A< 2) 

IF ( N , EQ i 2 ) G 0 TO 35 
DO 30 K s 3,N 

Y = T»Y«-AtK> 

32 CONTINUE 
35 P < I ) B Y 

4g CONTINUE 

RETURN 

52 OO 80 I=1,M 
T 3 X { I ) 

Y=A(N)*T*A<N-1) 

IF <N,EU,2) GO To 75 
no 7D K® 3 , N 

Y=Y»T+ A (N+l-K) 

7? CONTINUE 
75 PU>=Y 

82 CONTINUE 
RETURN 
ENO 





00.00 
00203 
00300 
00400 
00500 
00600 
00700 
00800 
00903 
01000 
01100 
01200 
01300 
01400 
01500 
01600 
01700 
01800 
01900 
02003 
02100 
02203 
02303 
02403 
02503 
02603 
02703 
02803 
02903 
03003 
03103 
03200 
03303 
03403 
03500 
03600 
03703 
03800 
03900 
04303 
04103 
04203 
04300 
04400 
04503 
0 4603 
04700 
04803 
04903 
05003 
05103 
05200 
05303 
05400 
0550Q 
05600 
05703 
05800 
059 0l ?l 
06003 


SUBROUTINE GAUSStA, R,N, ABINf ,TOL) 

C SUBROUTINE Fo« THE INTEGRAL FROM A TO B 8Y GaUSSJAN QUADRATURE 

c h c j ; l > are T H e heroes or t h e cl*i)ST lece n d r e polynomial 

c W(j;l> ARE The CORRESPONDING WEIGHTS 

c calls subroutine funct which defines the integrand 

C ABlNT IS THE OUTPUT 

implicit real *b <a-h#o* 2> 

DIMENSION 2(5,9)# W ( 5 # 9 ) # X ( 10 ) , F ( 10 > 

‘DATA ?/ .‘577350 26918962600, 4*0.00, 

1 0.03, .77459666924146303# 3*0,00# 

2 .33998134358485603, .86113631159405300, 3*0;D0, , 

3 0VO0, ,53846931013568303, .90617984593866403# 2*0,00. 

4 .23«6l9l86083l97D3, , 66 1209 3864 66 265q0 , , 932469514203 152D0 # 3* 

5 0703, .40584515137739700, ,74153118559939400# 

6 .94910791234275903, 3.00, 

7 .18343464249565003# , 525532 4 0991 632 9 0& # 1 796666477 41362700 # 

8 .96328985649753600# 0700, 

9 0VO0, .32425342340383900, , 61337143270059000 i 
A .83603110732663603# , 9681 60239507 626D0 # 

B .14887433098163103, . 4 33395394 12924 7 q 0 , ,67940956029902400# 

C .86506336668898503# ,97393652831717200/ 

DATA W/1,00, 4 *0 • P0 i 

1 . 88888 80888 9888903 #7 55555555555555600 1 3*0 .' D0# 

2 . 652145154862546D3, , 34785434513745400 , 3*0700# 

3 . 56080808988888903 # , 47862867049936600 # 

4 .236926885056189D3# 2*0 '.’D 3# 

5 .46791393457269103# ,36076157304813900# 

6 .17132449237917003, 2*0700# 

7 ,41795918367346903# ,36183005350511900, 

8 . 279705391469277 DC, ,1294849661688700^, 0700 i 

9 .36260378337836203, .31370664587788700, > 

A .222301034453374O0, .10122853629037600, 0700# 

B .3 30 23 9 3 5500 126 (’00, ,31234707704030300, .26361069640293500# 

C .18^64816369485700, , 08127438336157400# • 29552422471475300^ 

0 ,26926671930999600# ,21908636251598200# 

E , l4945l349i50581D0# , 066671344308688Q 0 / 

FACHs(B-A)/ 2 ,O 0 

IF (OARS(FACM) .GT.TOL) CO TO 5 
ABINT = 0.O 0 
WRITE (6,2) 

2 FORMAT <* (B^A) .LT.TOLI ABINT SET TO HERO » j 
RETURN 

5 F ACPs ( 9* A ) /2 »D0 

AB I NT *0 » 00 
EN = N 

EN02*EN/2 , 00 
W02=EN02 

!F( (N-2*N02) ,EQ,0) Kl»2 • 

!F((N-2*NO2),GT,0) Kl»i 

K2=ND2*1 

K 3 s 2 *K 2*1 

DO 10 K=Kl#K2 

t EpMsF AC M * i? ( K # N) 

X(K)sFACP-TERM . 

X ( K3 *.K) S FaCP^TERM 
xv continue 

IF (kl.EO.i) GO TO 15 
X( 1 )sEACP 
15 L r N ♦ 1 

CALL FUNCT<X,F,L> 



0610$ 


no 20 K=Kl,K2 

06203 


AB !NT3AB!NT+<nK)*F(K3*K ) >*W(K,N> 

06303 

22 

continue 

06403 


fF(KitEoVi> GO TO 25 

06500 


abintsabint^f<i)#w<i;n> 

06600 

25 

A8INtsABlNT#FACM 

06700 


RETURN 



